function [ output ] = inversion( data , Params , options  )

st   = data.st ;
s   = data.s ;
f   = data.f ;
W   = data.W ;
pop = data.pop ;
u   = data.u ;
area = data.area ;
   
% Construct steady-state unemployment rate
uss = ( s + Params.Delta ) ./ ( s + Params.Delta + f ) ;

% Recover zeta, v and bpv
run InvertZetaV.m

% House price contributions
I1 =   Params.omega .* u .* Params.b  ...
    .* W ./ bpv ./ ( 1 - Params.beta + Params.beta .* Ez ) ;
I2 = Params.omega .* ( 1 - u ) .* W ;
I3 = Params.psi .* W .* ( Ez ./ ( 1 - Params.beta + Params.beta .* Ez ) ) ;

% House prices
HousePrice = ( pop .* ( I1 + I2 + I3 ) ./ area ).^( 1 / ( 1 + Params.eta ) ) ;
    
% Productivities
L = W .* HousePrice.^Params.psi ...
    ./ ( bpv .* ( 1 - Params.beta + Params.beta * Ez ) )  ;

if strcmp(options.learning,'NoLearning') == 1
    lL = log(L) ;
    lH = zeros(size(L)) ;
elseif strcmp(options.learning,'Learning') == 1
    lH = log(    ( Params.mu + Params.Delta ) ...
              ./ ( Params.mu + Params.Delta + Params.phiLearning * uss ) ) ;
    lL = log(L) - lH ;
else
    error('Invalid options.learning')
end
lL  = lL - mean(lL) ;
lH = lH - mean(lH) ;

% Adjustment
J1 = Params.omega * u * Params.b ;
J2 = Params.omega * ( 1 - u ) .* bpv .* Ez ;
J3 = Params.psi * (1 - u ) .* bpv .* Ez ;
G = J1 + J2 + J3 ;

% Sorted zeta
zetaR = ( zeta - min(zeta) ) / ( max(zeta) - min(zeta) ) ;

% Mass of unemployed
Um = u .* ( 1 + u ) .* pop ;

% Assign output
output.zeta = zeta ;
output.v = v ;
output.bpv = bpv ;
output.Ez = Ez ;
output.HousePrice = HousePrice ;
output.lL = lL ;
output.lH = lH ;
output.FlowUNetAm = W ./ ( 1 - Params.beta + Params.beta * Ez ) ./ HousePrice.^Params.omega ;
output.G = G ;
output.zetaR = zetaR ;
output.Um = Um ;

end

